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Advance Earth System Science 
to meet the challenges of climate 
and environmental change. 
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NASA pioneered the interdisciplinary field of Earth System Science — the study of the 
Earth as an integrated system. This approach to studying the Earth as a single com- 
plex system is essential to understanding the causes and consequences of climate 
change and other global environmental concerns. Spaceborne instruments provide 
essential broad coverage, high spatial resolution, frequent sampling, and near-uniform 
accuracy and stability. Multiple on-orbit missions, including those flying in coordinated 
orbits as part of planned constellations, allow data to be acquired simultaneously 
on many important quantities, enabling investigations of the interactions among the 
coupled Earth processes that constitute the climate system. NASA’s research, cou- 
pled with that of our partners in the U.S. Global Change Research Program, provides 
much of the nation’s knowledge base for understanding, mitigating, and adapting to 
climate change. 
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• Definition 

“Exact pixel-to-pixel matching of two different 
images or matching of one image to a map ” 

• Multiple Source Data 

— Multimodal Registration 
— Temporal Registration 
— Viewpoint Registration 
— Template Registration 



What is Image Registration 
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Challenges in Image Registration 




• Remote Sensing vs. Medical or Other Imagery 

- Variety in the types of sensor data and the conditions of data acquisition 

- Size of the data 

- Lack of a known image model 

- Lack of well-distributed “fiducial points” resulting in lack of algorithms validation 


• Navigation Error 

• Atmospheric and Cloud Interactions 



Three Landsat images over Virginia acquired in August, October, and November 1999 
(Courtesy: Jeffrey Masek, NASA Goddard Space Flight Center) 




Atmospheric and 
Cloud Interactions 

Baja Peninsula, 
California; 4 different 
times of the day (GOES-8) 

(Reproduced from Le Moigne & 
Eastman, 2005) 
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Multitemporal 

Effects 

Mississippi and Ohio 
Rivers before & after 
Flood of Spring 2002 
(Terra/MODIS) 

(Reproduced from Le Moigne & 
Eastman, 2005) 




Relief Effect 

SAR and Landsat-TM 
Data of Lope Area, 
Gabon, Africa 

(Reproduced from Le 
Moigne & Eastman, 2005) 
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Image Registration or Precision Correction 



Navigation or Model-Based Systematic Correction 

- Orbital, Attitude, Platform/Sensor Geometric Relationship, Sensor 
Characteristics, Earth Model, etc. 


• Image Registration/Feature-Based Precision Correction 

— Navigation within a Few Pixels Accuracy 

— Image Registration Using Selected Features (or Control Points) to Refine 
Geo-Location Accuracy 

• Image Registration as a Post-Processing or as a Feedback to 
Navigation Model 



• (Towsnhend et al, 1992) and (Dai & Khorram, 1998): small error in registration 
may have a large impact on global change measurements accuracy 

• e.g., 1 pixel misregistration error => 50% error in Vegetation Index (NDVI) 
computation (using 250m MODIS data) 



Human-induced land cover changes observed by Landsat-5 in Bolivia in 1984 and 
1998(Courtesy: Compton J. Tucker (Hid the Landsat Project, NASA Goddard Space Flight 

Center) 

• Influence of image registration on products validation 

• Impact of misregistration on legal, economic and sociopolitical (e.g., resource 
management), etc. 
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Image Registration Frameworks 



Mathematical Framework 


- Il(x,y) and I2(x,y): images or image/map 

- find the mapping (f,g) which transforms II into 12: I2(x,y) = g(Il(fx(x,y),fy(x,y)) 

» f : spatial mapping 
» g: radiometric mapping 

- Spatial Transformations “f ’ 

- Translation, Rigid, Affine, Projective, Perspective, Polynomial, ... 

- Radiometric Transformations “g” (Resampling) 

- Nearest Neighbor, Bilinear, Cubic Convolution, ... 


• Algorithmic Framework (Brown, 1992) 

1. Feature Extraction 

2. Feature Matching 

3. Image Resampling 
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NASA Goddard Image Registration Group 





• 1 994: First results on the utilization 
of orthogonal Daubechies wavelets 
for image registration 


Figure 1 
Original Image 
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Figure 2 Figure 3 

Wavelet Coefficients Corresponding Wavelet Coefficients Correspond 

to Figure 1 to Figure 1 rotated 44 degrees 
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• Study of rotation- and translation-invariant wavelet 
filters (Spline, Simoncelli) 

• Study of different matching strategies and metrics 

• Parallel implementations (SIMD/MasPar, Beowulf 
Cluster, MIMD/Cray-T3E, FPGA-Hybrid) 

• Development of 
image registration 
framework based 
on Brown’s 
framework 




NASA Goddard Image Registration Group 
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Experiments ... Datasets (1) 






• Synthetic Data Experiments 
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Experiments (1) ... Analysis Samples 



Various Features; Convergence as a function of noise 
and radiometric variations 


(white areas - regions of convergence with errors less than threshold, e.g. 0.5) 


tAU-SfflC TRU-SimB TRU-SifnL 



Warping + Noise 


1RU-$0C TRu-&n8 Tfiy-Siml, 



Warping + PSF 


• Simoncelli-based methods outperform 
Spline pyramid-based methods 

• Optimization based on Mutual Information 
does not perfom better than L2-Norm 

• Simoncelli-LowPass better than Simoncelli- 
BandPass for Low Noise and Same 
Radiometry and for Initial Guess Sensitivity 


TRO-SpC TRU-SwnB TRU-SJnK. 



Warping + PSF + Noise 







Experiments ... Datasets (2) 



• Multi-Temporal Data 


- Landsat-5 and -7 (chips and corresponding windows) 




1 Landsat chip 

and 4 corresponding windows 


7 Landsat chips 



• Multi-Sensor Data 


- EOS Validation Core Sites 

- IKONOS/Landsat-7/MODIS/SeaWiFS 

• Red and NIR bands for each sensor 

• Spatial resolutions: IKONOS: 4m; ETM+: 30m; MODIS: 500m; 
SeaWiFS: 1000m 

- 4 different sites: 

• Coastal Area: VA, Coast Reserve Area, October 200 1 

• Agriculture Area: Konza Prairie in State of Kansas, July to 
August 200 1 

• Mountainous Area: Cascades Site, September 2000 

• Urban Area: USDA Site, Greenbelt, MD, May 2001 
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Experiments ... Datasets (3) 



Multi-Sensor Data 



ETM/IKONOS - Coastal 

VAData 



ETM/IKONOS - 
Agricultural Konza Data 
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Experiments (2 and 3) ... Analysis Samples 



Global transformation vs. manual registration (or “ ground truth ”) 
parameters for 4 Scenes in DC mutitemporal dataset 
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Self-Consistency Study of the Mutual 
Information Results 
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Number of cases that converge (out of 32) for the DC dataset, running 4 


algorithms and different features with the initial guess varying between the 
origin (d=0.0) and ground truth (4=1.0) 


GOAL: DEFINE A “REGION OF CONVERGENCE” AND A “REGION 
OF DIVERGENCE” FOR EACH ALGORITHM 

RECOMMENDATION FOR UTILIZATION OF ALGORITHMS AND 
ITS COMPONENTS 
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Web-based Image Registration 


Tutorial 
Contact Us 

1. Open Image files 

Input Image: ( Browse ) 

Reference Image: I ( Browse ) 

2. Select Registration Algorithm 

O ureg 

© TRUMI 
o SPSA 

O Use All Algorithms 

3. Customize Algorithm 

Select Parameters 1 j ) 

4. Perform Registration 

( Register Images ") f Reset ^ 
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Contact Us 


Successful Completion! 


Reference Image 



Input Image 



Parameters Selected 

Rows = 256; Columns = 256; Wavelet Type = Spline 


Result: 


Output Image 


Jransletion coefficients: 


floatation angle = 4.000000 

Shift.x = 5.0000CC 1 

Shifty = 2.000000 

Scale.* = 1.000000 1 

1 jgfrjg) 

Scale.* = 1.000000 

f i ^ jHBh 


• Save Output Image ' ' Email Result^ ' Register This Image Set Again N \ Register Another Image Set ''l 


Copyright 2005. Code 606. 3. NASA/Goddard Space Flight Center. Greenbelt. Mary land, USA. 
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• Image Registration for Remote Sensing , ed. J. Le Moigne, N.S. Netanyahu and 
R.D. Eastman, Cambridge, UK:Cambridge University Press 


• Foreword by Jon A. Benediktsson 

• Contributors: S. Baillarin/CNES; D.G. Baldwin/Univ. of Colorado; M. 
Bemard/SPOT Image; A. Bouillon/Institut Geographique National; J.L. Carr/Carr 
Astronautics; R. Chellappa/UMD; Q-S. Chen/Hickman Cancer Center; A. Cole- 
Rhodes/Morgan State Univ.; R.I. Crocker/Univ. of Colorado; R. Davies/Univ. of 
Auckland; D .J. Diner/NASA JPL; W.J. Emery/Univ. of Colorado; A.A. 
Goshtasby/Wright State Univ.; V.M. Govindu/Indian Institute of Science; V.M. 
Jovanovic/NASA JPL; C.S. Kenney/UC Santa Barbara; B.S. Manjunath/UC Santa 
Barbara; J. Morisette/USGS; D.M. Mount/UMD; M. Nishihama/Raytheon 
@NASA GSFC; F.S. Patt/SAIC @NASA GSFC; S. Ratanasanya/form. UMD; K. 
Solanki/UC Santa Barbara; H.S. Stone/form. NEC Research Lab; J. Storey/SGT 
@USGS; S. Sylvander/CNES; B. Tan/ERT @NASA GSFC; P.K. 
Varshney/Syracuse Univ.; R.E. Wolfe/NASA GSFC; C. Woodcock/Boston Univ.; 
M. Xu/Syracuse Univ.; I. Zavorin/form. UMBC@NASA GSFC; M. Zuliani/UC 
Santa Barbara 



• Part 1 - The Importance of Image Registration for Remote 
Sensing 


• Part II - Similarity Metrics for Image Registration 


• Part III - Feature Matching and Strategies for Image 
Registration 


• Part IV - Applications and Operational Systems 


Part V— Conclusion and the Future of Image Registration 
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F eature Matching 

Feature (Extraction), Similarity Metrics, 
Transformations, and Matching Strategies 


Nathan S. Netanyahu 
Dept, of CS, Bar-Ilan University, Israel, 
and CfAR/UMIACS, Univ. of Maryland 
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• Given a reference image, Ifx, y), and a sensed image fix, y), find 
the mapping (T , g ) which “best” transforms f into I 2 , i.e., 

I 2 (*> y) = SifiT p (x, y), T p (x, y))), 

where T denotes spatial mapping and g denotes radiometric 
mapping. 

• Spatial transformations: 

— Translation, rigid, affine, projective, perspective, polynomial 

• Radiometric transformations (resampling): 

— Nearest neighbor, bilinear, cubic convolution, spline 
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• Gray levels 

• Salient points 

— Edge-like, wavelet coefficients (Simoncelli and 
Freeman ‘95) 

- Comers (Keamy et al. ‘87, Harris and Stephens ’88, 

Shi and Tomasi ‘ 94) 

• Lines 

• Contours, regions (Govindu et al. ‘99) 

• Scale invariant feature transform (SIFT), Lowe ‘04 



Feature (Extraction) 
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— Minimize the sum of squared errors (SSD) over overlapping 
subimage 

M-lN-1 

SSD(x, y ) = SB' 1 (m, ri) — I 2 (m — x,n — y)f 

m = 0 n=0 
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Similarity Metrics (cont’d) 


Cross-correlation 

— Maximize cross-correlation over image overlap 

M - 1 N - 1 

/, (x, y) o I 2 (x, y) = (x + m,y + n) 

m = 0 «=0 



Normalized cross-correlation (NCC) 

— Maximize normalized cross-correlation 


M - 1 JV-1 



NCC h ,,S x ^y) = 


/, (m, n) — /[ 1 f/ 2 (x + m, _y + n) — 1 2 1 
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— -|2 
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— -|2 
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m = 0 / 7=0 


772—0 72—0 
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Similarity Metrics (cont’d) 




• Mutual information (MI): 

Maximizes the degree of statistical dependence between the images 


Mi{i v i 2 )=YT.Pk,iM’sM°% 


g\ § 2 


V 


P/,./ 2 (gpg 2) 

P U (g, (g 2 ) 
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or using histograms, maximizes 


m/(/ 1? / 2 ) 


\_ 

M 


ZIXJg.^Hog 

g\ g 2 


f (gpgz) ^ 

V ^,(gl)-\(g2) y 


where Mis the sum of all histogram entries, i.e., number of 
pixels (in overlapping subimage) 
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• Partial Hausdorff distance (PHD): 


H K (l\>h) = K'L, dist ( A , p 2 ) , 


where 1 < K < I/, I (Huttenlocher et al. ‘93, Mount et 





Similarity Metrics (cont’d) 


Discrete Gaussian mismatch (DGM) 


w (a ) = exp 


dist(a,/ 2 ) 2 ^ 




2(7 
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where w (a) denotes the weight of point a, and 


DGM (7 


M = i-^TTT — 

I A I 


is similarity measure ranging between 0 and 1 (Mount et al., Ch. 8) 


Transformation Functions 


Translation-only, rigid 

Rotation, scale, and translation (RST) 

Affine (6 degrees of freedom) 



x' = scos^-x-ssin#- + 
y ' = s sin 6 ■ x + s cos 0-y + t 



r s cos 6 
ssin# 




0 


-s sin 6 tS\ 

A 

s cos 6 t 

0 lj 


Projective/homography (e.g., for perspective effects in image 
mosaicing; Govindu and Chellappa, Ch. 10); 8 parameters 
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• Weighted linear transformation (Goshtasby, Ch. 7); adaptive 
transformation, continuous and smooth, applied to multiview 
images with local geometric differences, and maps an entire 
image to another 

— Interpolating surface is a weighted sum of planar patches, 
each of which passes through a control point and provides a 
desired gradient, i.e., 



Transformation Functions (cont’d) 


f(x,y) = 


ZwA(*,y)A(*>y) 


Zw 

for monotonically decreasing weight R j (x,y) = [(.v-.v f + ( r — r r] 


and 


A y) = a , ( x ~ x i ) + b, (y - y, ) + F t 
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Transformation Functions (cont’d) 




Reference Sensed 


Registered 



Source: Goshtasby, IR Tutorial, CVPR 77 
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IR Components (Revisited) 



Features 


Similarity 

measure 


Matching 

strategy 


Gray levels 


Spline or Simoncelli 
LPF 


Simoncelli 

BPF 



Thevenaz, 

Ruttimann, 

Unser 

optimization 


Spall’s 

optimization 
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• Exhaustive search (exponential in dimensionality of space) 

• Fast Fourier transform (FFT) 

• Numerical optimization (e.g., steepest gradient descent wrt SSD, 
NCC, and MI (Thevenaz, Ruttimann, and Unser (TRU) ‘98; Spall 
‘92)) 

• Robust transformation estimate (e.g., RANSAC, LMS) if (most) 
correspondences are known (via SIFT-like) 

• “Correspondenceless”, e.g., correlation of descriptor 
distribution/feature consensus (Govindu et al. ‘99) 

• Robust feature matching (RFM), e.g., efficient subdivision and 
pruning of transformation space; Huttenlocher et al. ‘93, Mount et 
al. ’99, Netanyahu et al. ‘04 



Matching Strategies 
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• Frequency domain-based approach (Stone, Ch. 4) 


— Efficient computation of correlation as inverse of 
F*(u,v)F 2 (u,v) 


- Practical implementation (extension to NCC, masking 
invalid pixels, optimized computation) 

- Finding (small) rotational and scale differences (by 
matching chips) 

- Subpixel registration for translation-only using phase 
estimate (also in case of image aliasing) 

- Rotation and scale estimate by casting to log-polar 
coordinates 
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Matching Strategies (cont’d) 


• Matched filtering (Q. Chen, Ch. 5) 

— Maximize SNR (using theory of linear systems) 

— Apply phase-only and symmetric phase-only matched filters for 
translation-only IR 



Phase product = 


F*(u,v) F 2 (u,v ) _ 

_ g-jiMtjc+Vty) 

F t (u,v) F 2 (u,v) 


- Apply Fourier-Mellin transform for rotation and scale changes; 
transform represents these parameters as translational shifts in 
log-polar coordinates of magnitude of Fourier spectrum, i.e., 
first estimate rotation and scale, followed by translation estimate 




Matching Strategies (cont’d) 



Rotation and 
scale estimate 



Translation 

estimate 
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Pair of SPOT images and their registration, using symmetric 
phase-only matched filters on their Fourier-Mellin transforms 
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Matching Strategies (cont’d) 



Numerical optimization (Cole-Rhodes and Varshney, Ch. 6; 
Cole-Rhodes and Eastman, Ch. 12) 

— Powel’s, Brent’s (1-D), simplex, etc. 


- Steepest descent/ascent variants 

•Standard P*+i = Pi ~ K%k 

• Newton-Raphson P*+i = P * - \ H l §/. 

• Levenberg-Marquardt p t+1 = p k -(K k + ^diag[H t ])"' g k 

— Apply to various similarity metrics, e.g., SSD (Eastman and Le 
Moigne ‘01), MI, etc. 

» Explicit computation of gradient (and Jacobian/Hessian), e.g., 
Thevenaz and Unser ‘00 

» Stochastic approx. (Spall ‘92); Cole-Rhodes etal. ’03; Cole-Rhodes 
and Varshney, Ch. 6 
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Matching Strategies (cont’d) 




Pair of Landsat images over DC 



MI surfaces of above (level 1 and 4) images, using B-spline interpolation 
(Cole-Rhodes and Varshney, Ch. 6) 
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Matching Strategy (cont’d) 




• Alignment via local geometric distributions 

(Govindu and Chellappa, Ch. 10) 




-150 -100 -50 0 50 1 00 150 -150 -100 -50 0 50 100 150 

(c) 



Slope angle distributions and their correlation 
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• Robust feature matching (RFM) (Mount et al Ch. 8) 

- Space of affine transformations: 6-D space 

- Subdivide-. Quadtree or kd-tree. Each cell T represents a set 
of transformations; T is active if it may contain t ; o/w, it is 

killed 

— Uncertainty regions (UR’s): Rectangular approximation to 
the possible images r(a) for all r e T,a e /, 

— Bounds: Compute upper bound (on optimum similarity) by 
sampling a transformation and lower bound by computing 
nearest neighbors to each UR 

- Prune : If lower bound exceeds best upper bound, then kill the 
cell; o/w, split it 



Matching Strategy (cont’d) 



RFM-based registration of Landsat images over DC using wavelet 
features and PHD similarity measure (Netanyahu et al. 4 04) 
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• Computational efficiency 

— “Culling” feature points via, e.g., condition theory 
(Kenney et al. ‘03, Ch. 9) 

— Efficient numerical or discrete algorithmic 
procedures 

— Hierarchical pyramid-like (wavelet) decomposition 

— Use landmark chip database (instead of a large 
scene) or alternatively, extract automatically 
corresponding regions of interest using 
mathematical morphology (Plaza et al. ‘05, ‘07) 



Matching Strategy (cont’d) 
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• Use Cramer-Rao bounds as performance benchmark 
for performance evaluation of image registration (Xu 
and Varshney, Ch. 13) 



Miscellaneous 



From Theory to Practice 
Operational Requirements 

Roger D. Eastman 

Loyola University, Baltimore, 

Maryland 
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Why isn’t this problem solved by now? 



• A wealth of approaches ! 

• SIFT, ASIFT, BSIFT, SIFT/NCC, SIFT/FLOUR 


• Beat the problem to death with terminology 

• “Assume we have a Banach space ...” 

• Many smart people wielding heavy mathematical weapons 
against a relatively fixed problem- why hasn’t the problem 
yielded? Why no gold standard algorithm? 



Operational Satellite Team s solve it every day 


•GOES -Carr, Chapter 1 5 
•MISR - Jovanovic et al, Chapter 1 6 
• AVHRR - Emery et al, Chapter 1 7 
•Landsat, Storey, Chapter 1 8 
•SPOT, Ballarin, Chapter 19 
•VEGETATION, Sylvander, Chapter 20 
•MODIS, Wolfe et al, Chapter 21 
•SeaWiFS, Patt, Chapter 22 


And it’s often solved the old-fashioned way (200 
___ ___ Normalized C 


Instrument 

Satellite 

Resolution 

Similarity 

Subpixel 

ASTER 

Terra 

15m-90m 

NCC w/ DEM 

Fit to surface 

GOES 

GOES l-M 

1 km-8km 

NCC w/ vector 
coastlines 

Bi-section 

search 

MISR 

Terra 

275m 

NCC w/ DEM 

Least squares 

MODIS 

Terra 

250m-1km 

NCC w/ DEM 

Fixed grid 

HRS 

SPOT 

2.5m 

NCC w/ DEM 

Not described 

ETM+ 

Landsat-7 

15m-60m 

NCC to arid 
region CPs 

Fit to surface 

VEGETATION 

SPOT 

1 km 

NCC w/ DEM 

Not described 



• Geodetic accuracy 

- Database of GCPs derived from USGS data 
— Normalized correlation 

- Updates navigation models 
— Results: RMSE ~54m 

• Band-to-band registration 

— Selected tie-points in high-freq. arid regions 
— Normalized correlation 

- Subpixel by second order fit to 3x3 neighborhood 

- Result: 0.1 to 0.2 subpixel 
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• Know models of sensor/platform/ 

• Have access to complete data set 

• Have continuing demands/responsibility 

• Are registering same plots of land again and 
again - can invest effort in data preparation 

• Can’t take big risks on unproven methods 



Operational teams requirements 
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Know platform: Landsat team knowledge 


Sensor geometry 

— Band to band 

Sensor to platform 



— Sensor to sensor 

• Orbit 

— Platform to Earth 

• Terrain data 

-DEM 

• Radiometric model 



+X 

tffifctrations USGS/NASA 



Image Registration 
for Remote Sensing 


Invest in data: 


ETM+ Chip 




Image Registration 
for Remote Sensing 


Know data: GOES channel 1 (Baja) 



• Contrast reversal 
day to night 

• Requires use of 
contour matching 
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Use DEMs: Digital Terrain Models 



Taking terrain into account in matching 
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Use proven methods: Landsat 7 library 


Clean data, go fast 

— Use Normalized Grey-Scale Correlation 

Missing data/gaps, need robustness 

— Use Mutual Information 



• Available alternative 

— Use Robust Feature Matching 



• Know what data is for 

• Have to fuse many data sets 

• Have access to ancillary data 

• Know cultural and historical data 

• Don’t need one magic method - need toolbox 
of many approaches 
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• Different communities/literature/requirements 

— Photogrammetry 
— Computer vision/image processing 
— Operational teams 

— Remote sensing/Earth scientists/end users 

• Demanding/varying mission requirements 

— Caution in system design, new methods 

• Expensive sensors and images 

— Hard to share data or complete models 



Institutional challenges to solving IR for RS 



Conclusion 


Jacqueline Le Moigne 
Nathan S. Netanyahu 
Roger D. Eastman 



• Satellite sensing/imaging in full expansion 

— Explosion of commercial satellites 

— Exploring distant planets (Moon, Mars, etc.), e.g. Lunar Reconnaissance Orbiter 
(LRO) 

• Future research and challenges 

— Combining multiple band-to-band registrations (e.g., hyperspectral data) 

— Automatically extracting windows of interest (decreasing processing time and 
increasing accuracy) 

— Dealing with other data sources (e.g., planetary imagery, or verification of 
optical systems) 

— Integration and fusion of multiple source imagery (various satellites, vector map, 
airborne, ground data, etc.) 

— Onboard implementations on specialized hardware 

— Multistage registration algorithms combining multiple principles and approaches 
and utilizing interdisciplinary systems engineering approach , thus increasing 
algorithms robustness and applicability 
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The Autonomous Land Vehicle (ALV) Project 





Thank You! 


Jacqueline Le Moigne 
Nathan S. Netanyahu 
Roger D. Eastman 


